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SUMMARY 


A numerically generated orthogonal coordinate system (with the body surface 
and shock wave as opposite boundaries) has been applied with a time asymptotic 
method to obtain steady-flow solutions for axisymmetric inviscid flow over 
several blunt bodies including spheres, paraboloids, ellipsoids, hyperboloids, 
hemisphere cylinders, spherically blunted cones, and a body with a concavity 
in the stagnation region. Comparisons with experimental data and the results 
of other computational methods have demonstrated that accurate solutions can 
be obtained with this approach. The numerically generated orthogonal coordinate 
system used in the present paper should prove useful for applications to complex 
body shapes, particularly those with concave regions. In addition, the use of 
the present orthogonal coordinate system simplifies the form of the governing 
equations and simplifies the application of boundary conditions at the body sur- 
face and shock wave. 


INTRODUCTION 

Previous investigations of the direct blunt-body problem have demonstrated 
the utility of treating the bow shock wave as a discrete discontinuity. (See, 
for example, refs. 1 to 5.) In each of these studies, some form of coordinate 
transformation was used to map the physical space within the shock layer into a 
rectangular computational domain, with the bow shock wave as one boundary and 
the body surface as the opposite boundary. These transformations generally pro- 
duce a nonorthogonal coordinate system in the physical space which complicates 
the form of the equations of motion. However, they usually produce a uniform 
grid system for calculating the flow-field solution and simplify the application 
of the boundary conditions at the shock wave and the body surface. 

Recently, Graves (ref. 6) applied a method for numerically generating an 
orthogonal coordinate system between two arbitrary continuous curves. By allow- 
ing one of the curves to be the body surface and the other to be the bow shock 
wave (which can move with time) , this method can be used to construct an orthog- 
onal coordinate for solving the blunt-body problem. An orthogonal coordinate 
system of this type leads to a much simpler form of the equations of motion 
and further simplifies the application of boundary conditions at the body sur- 
face and the shock wave. In addition, this type of coordinate system can be 
applied to complex body shapes, including bodies with concave regions. 

The purpose of the present paper is to explore the application of the 
numerically generated orthogonal coordinate system presented in reference 6 to 
flow-field problems. For this purpose, the inviscid, axisymmetric supersonic 
flow over a blunt-nose body will be investigated. 



SYMBOLS 
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B b 
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d 

h-| ,h2 

h 3 

I 
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A 

L 

M 

P 

*b 

R 

s 

t 

u 

V s 


V n 

V t 


^ /* 2 

speed of sound, a/V^ ; also length of Z-axis for ellipse 

bluntness parameter 

length of r-axis for ellipse 

specific heat at constant pressure 

specific heat at constant volume 

✓V f<\ 

diameter of cylinder, d/L 

metric coefficients in and redirections, h-j/L and h 2 /L, 
respectively 

xn y /\ 

metric coefficient in circumferential direction, h 3 /L 

^ ^ 2 

static enthalpy, I/V^ 
index for ^-direction, E, = (i - 1) 
index in redirection, r| = ( j — 1 ) Ari 
dimensional characteristic length, m 
Mach number 

/v A ~ 2 

pressure, P/P^ 

-A j s\ 

body nose radius, R^/l 
radius defined in figure 1, r/l 
coordinate arc length, s/L 
time, Ev^/L 

velocity component in ^-direction, u/v^ 

A i A 

shock velocity, Vg/V^ 
free-stream velocity, v/v^ = 1 
dimensional free-stream velocity, m/sec 
component of velocity normal to shock wave, V n j 

/\ I A 

component of velocity tangent to shock wave, Vt/V oo 
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V 


/ A 

V co 

Z,r Cartesian coordinate system (see fig. 1), z/£, ?/l 

Zp location of pole for coordinate (see fig. 1), Z p ^£ 

8 shock wave angle 

Y ratio of specific heats c p/ /c v 

AR shock-layer thickness in polar coordinates, Ar/£ 

At time increment 

An increment in n-coordinate 

AEj increment in ^-coordinate 

C,n orthogonal coordinate system (see fig. 2) 

p density, (3/p^ 

<j> angular coordinate defined in figure 1 

Superscripts: 

(~) dimensional quantity 

(”) predicted quantity 

Subscripts: 

max maximum 

o stagnation 

s shock 

00 free stream 


METHOD OF SOLUTION 

The problem is to compute the inviscid, axisymmetric supersonic flow over 
a blunt- nose body. Solutions are to be obtained throughout the subsonic region 
and downstream into the supersonic region (fig. 1). 
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Coordinate System 


Using the method described in reference 6, a general orthogonal coordinate 
system was numerically constructed between the body surface (n = 0) and the bow 
shock wave (U = 1) extending from the stagnation streamline (£ =0) to a down- 
stream outflow boundary (£ = 1) located within the supersonic portion of the 
shock layer (fig. 2) . The line element for this orthogonal coordinate system 
(£,U) is given by the equation 


ds^ = hi 2 d£2 + h2 2 du^ 


(1) 


where hi is the metric coefficient in the ^-direction and h 2 is the metric 
coefficient in the ri-direction. 

The ^-coordinates can be chosen so that grid points are equally spaced 
along the body surface in the physical plane or so that they may be concentrated 
in regions of large gradients. Similarly, the n-coordinates can be chosen so 
that grid points are equally spaced between the body and the shock wave or so 
that they may be concentrated in any desired region. However, in both cases 
the grid system is uniform in both the £- and n-directions in the computational 
plane (fig. 3) . Any nonuniformity in the physical plane is taken care of by the 
metric coefficients hi and h 2 , which are calculated each time the coordinate 
system is calculated. 


Flow-Field Equations 

The partial differential equations governing the axisymmetr ic compressi- 
ble flow of an inviscid nonconducting fluid in general orthogonal coordinates 
are 

Continuity : 

1 3u 

h l 9? 

v / u 9 ^3 v 

h-j h2 3ri h3 3£ ^2^3 


9p u 3p v 9 P 

— + — — + — — + p 

9t hi 9£ h2 9f| 


1 3v u 9h 2 

h 2 9u h 1 h 2 9 £ 


3hs^ \ 
3rt j 


= 0 


( 2 ) 


4 



^-momentum ; 


3u u 3u v 3u / v ^ h 2 u 1 /I 3p\ 

3t + h-, 3£ + h 2 3n V \h 1 h 2 3? h 1 h 2 3n / P\h 2 3?/ ” 


(3) 


ri-momentum : 


3v u 3v v 3v / v u ^m\ 1 /l 3p\ 

— + — — + — — + ul - + — ( — — ) = 0 

3t h-| 3£ h 2 3ri yM^ 2 ^ h-|h 2 ^ / p\h 2 


(4) 


The preceding system of equations in general orthogonal coordinates is 
much simpler than the system of equations that is obtained if the space between 
the body and the shock wave is mapped into a rectangular computational domain 
(refs. 1 to 3), where the resulting coordinate system in physical space is 
nonorthogonal. 

In order to close this system of equations, the integrated form of the 
steady-flow energy equation 


1 co + 


u 2 + v 2 


= I, 


(5) 


is used along with the ideal gas equation of state 


P = 



( 6 ) 


Steady-flow solutions are obtained by integrating the system of equa- 
tions (2) to (6) in time, from an assumed set of initial conditions, until 
steady state is reached. Since the integrated form of the energy equation 
(eq. (5)) is being used, the results are not accurate during the transient 
phase of the solution. However, equation (5) becomes more accurate as the 


5 


steady state is approached and the correct converged steady-flow solution is 
obtained. Using the energy equation in this form speeds up the computational 
procedure. 

Along the stagnation streamlines u = h 3 = 0 and 3h 3 /3£; = h-j ; thus the 
u 9 **3 v 9 ^3 

terms and appearing in equation ( 2 ) are of the indefinite 

h-j I 13 3£ ^ 2^3 9 *1 

form (0/0). Taking the limit of these terms as £ -*■ 0 by using L'Hopital's 
rule 


( u 
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&-0 \ h 1 h 3 
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( 8 ) 


Thus on the stagnation line (£ = 0) , the limiting form of the continuity equa- 
tion ( 2 ) becomes 


3p v 3p 1 

— + + p — 

3t h 3 3r| h 3 


3v A 3u v 9h 1 
3ri yh-j 3£ h-j h 3 t 


(5 = 0) (9) 
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Equations (2) to (9) are written in nondimens ional form using the following 
definition for the nondimensional variables: 


s\ 

U /Voo 

~ A 

= ty L 


h l = h l/ L 

P/Pco 

/\ f * N 

h 2 = h 2/ L 


h 3 = h 3 /£ 

p/pjj 



Finite-Difference Scheme 

The solution of the flow field is obtained by integrating the compressible, 
time-dependent Euler equations described in the previous section. The integra- 
tion is carried out using the explicit finite-difference scheme of Brailovskaya 
(ref. 7) . 

Using the following computational module 



this two-step difference scheme can be represented as follows: 
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Predictor step : 


Ola) 


01b) 


where F can represent p, u, or v. The time derivatives appearing in equa- 
tions (11a) and (lib) are obtained from equations (2) to (4). The spatial deriv- 
atives in equations (2) to (4) are evaluated using central difference formulas 
except at the boundaries where three-point, one-sided difference formulas are 
used. 


t+At ♦. / 9F Y- 

p i.j - p lj - b 


VVi.j 


Corrector step : 


t+At 

r i#j 


= F 




/ 9 F \t+At 

+ At UA,j 


The numerical stability limit for this difference scheme is the well-known 
CFL condition given by the equation (ref. 8) 


(At) 


max 


M 


+ a, 


h-t AC h 2 An y (h-| AC ) 2 (h-! An ) 2 


( 12 ) 


The local maximum time step at each grid point is calculated by using equa- 
tion (12), and then a fraction of the local value (usually 50 to 80 percent) 
is applied at each grid point to advance the solution in time. This leads to 
an inaccurate representation of the time-dependent nature of the solution, but 
it has been found to accelerate convergence by a factor of approximately 1.5 
to 2 when compared with using the minimum global time step. 


Boundary and Initial Conditions 

The boundary condition of no flow normal to the body surface is 


v = 0 


(n = o) 


(13) 


since the body surface is the inner boundary of an orthogonal coordinate system. 
The remaining properties at the surface are computed using three-point forward 
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differences to calculate the spatial derivatives in both the predictor and cor- 
rector steps of the finite-difference scheme described in the previous section. 

At the outflow boundary (£ = 1) and the shock wave (n = 1), three-point 
backward differences (away from the boundary) are used. At the shock wave, the 
pressure thus calculated is used with the shock-wave relations (described in 
the next section) to calculate the shock velocity and then to update all other 
properties at the shock wave. 

Along the stagnation streamline (£ 


u = 0 

3v 3p 9p ^2 

Si ~ ii " 5i ~ aF 

and the remaining derivatives are computed by taking into account the asymmetry 
in the u-component of velocity. 

To start the solution, a shock-wave shape is assumed and the shock-wave 
velocities are set equal to zero. Properties at the shock wave are computed 
from the shock-wave relations. Properties on the body surface are computed 
from modified Newtonian theory, assuming that the entropy is constant on the 
surface. Properties at a constant value of E, across the shock layer are 
approximated by linearly interpolating between the values at the shock wave 
and the body surface. 


= 0) , the following conditions are set: 


04) 


Solution for Shock Wave 

A method similar to that used by Weilmuenster and Hamilton (ref. 9) has 
been used to track the movement of the shock wave during the transient portion 
of the solution. The method can be summarized as follows (details are given 
in the appendix) . First the pressure is computed on the downstream side of 
the shock wave by using the two-step difference scheme described previously. 
With the pressure at the shock wave p s known, the density is computed from 
the equation 


I eo + " (P S " PJ 


p s / Y \ 1 

M— ) ■ ; (Ps - p -> 


( 15 ) 
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and then the shock velocity is computed with the equation 



( 16 ) 


where V s and (V n ) are defined in the sketches in the appendix. Next, 
the movement of the shock wave is computed from the following differential 
equations: 


dz s 

— = V s sin B s (17a) 


dr s 

— = -V s cos 6 S (17b) 


where 

B s = tan -1 (dr s /dZ) 

by using the following two-step finite-difference approximations: 

Predictor step : 

(18a) 
(18b) 


(Zg) t+At = Zg + At(Vg SLU Bg) fc 
(r s ) t+At = r| - At(V s cos 3 S ) fc 
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Corrector step : 


(Z s ) t+At = z| + — k sin Bg) fc + (Vg sin B s ) t+At ] 

(r s ) t+At = r| f(V s cos B s ) fc + (V s cos B s ) t+At ] 

2 u J 


(19a) 


(19b) 


The coordinate system generator is used to update the shock wave after 
the predictor step and the entire coordinate system after the corrector step. 
As the solution approaches convergence, the shock velocity approaches zero 
and the shock wave remains essentially fixed in space. The shock velocities 
at convergence are generally less than 10~3 i n magnitude. 


Convergence and Computational Time 

The average change in density from one time step to the next was used to 
determine convergence. Based on the results of several test cases, it was found 
that the solution was converged when the average change in density was less than 
2 x 10“5. For the cases presented in the present paper, this generally required 
from 400 to 600 time steps. 

From timing studies on the Control Data CYBER 175 computer system, it has 
been found that approximately 3.4 x 10“^ sec/point/step were required to compute 
the flow field and approximately 7.8 x 10“^ sec/point/step were required to com- 
pute the coordinate system. Thus, in the present nonoptimized computer code, 
approximately 70 percent of the time is spent computing the coordinate system. 
One obvious method of improving the overall efficiency is to update the coordi- 
nate system less frequently. For example, it is possible to update the shock 
wave after each time step and hold the remainder of the coordinate system fixed 
for several time steps before recomputing. This would greatly reduce the per- 
centage of computational time spent on generating the coordinate system. 


Body Geometry 

Most of the body geometries that are considered in the present paper 
have a projection in the plane that is described by the equation for 

a general conic section 


r 2 = 2R b Z - B b Z 2 


( 20 ) 
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/N A 

where R b is the radius of curvature at z = 0 and B b is a bluntness param- 
eter. The bluntness parameter B b characterizes the eccentricity of the conic 
section. Its significance is better understood if it is noted that B b < 0 
generates a hyperbola, B b = 0 generates a parabola, and B b > 0 generates an 
ellipse (with B b = 1 for the special case of a circle) . For an ellipse, B b 
is related to the ratio of major to minor axis (b/a) by the relation 


B b = (b/a) 2 


( 21 ) 


Some of the body shapes considered in the present paper were spherically blunted 
cones. These body shapes were generated using equation (20) for the nose, with 
B b = 1 , followed by a straight-line segment tangent to the nose at the appro- 
priate cone angle. 

One body shape considered had a reverse curvature at the nose matched 
to a 30° conical afterbody at the point 5=0 and 9 = 1 and defined by 
the equation 


Z = - /3 r 2 + y5 r 3 


( 22 ) 


RESULTS AND DISCUSSION 

In this section, results of computation by the present method are compared 
with experimental data and with the results of other computational methods to 
demonstrate the applicability of the numerically generated orthogonal coordi- 
nate system of reference 6 for flow-field calculations. 


Sphere 

Pressure distributions over a sphere for y = 1.4 are presented in fig- 
ure 4 for free-stream Mach numbers of 2, 3, 4, 6.05, and 8.06. Experimental 
data from reference 10 and computed results from reference 3 are presented for 
comparison. The surface pressures calculated by the present method agree well 
with the other results for each Mach number considered. 

The shock-layer thickness for these conditions is presented in figure 5. 

The agreement with the other results is excellent except at M^ = 2 where the 
experimental data fall slightly below the calculated results. Since shock-layer 
thickness changes rapidly at low Mach numbers, this behavior could easily be due 
to a small change in the effective Mach number for the experimental tests which 
would change the shock-layer thickness but would have little effect on the mea- 
sured pressure distribution. 
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The calculated shock-wave shape and sonic-line location about a sphere at 
Y = 1.4 are presented in figure 6 for free-stream Mach numbers of 6.8 and 4.76 
Experimental data from references 11 and 12 are presented for comparison. For 
both conditions , the agreement with the experimental data is very good. 


Paraboloid 

The calculated surface pressure distributions over a paraboloid at 
Y = 1.4 are presented in figure 7 for Mach numbers of 3 and 10. Calculated 
results from references 3 and 13 are included for comparison. The results 
of the present calculation are in very good agreement with the other calcu- 
lated results. 

The shock-layer thicknesses for these conditions are compared in figure 8 
where, again, excellent agreement with the results of references 3 and 13 will 
be noted. 


Ellipsoid 

Pressure distributions for an ellipsoid with b/a = 0.5 and y = 1-4 are 
presented in figure 9 for Mach numbers of 4, 6.05, and 8.06. Experimental data 
from reference 10 and calculated results from reference 3 are presented for com 
parison. The present calculated results agree closely with the other data for 
each Mach number . 

Pressure distributions for an ellipsoid with b/a =1.5 and Y = 1 *4 are 
presented in figure 10 for Mach numbers of 3, 4, 6.05, and 8.06. Experimental 
data from reference 10 and calculated results from reference 3 are presented 
for comparison. The present calculated results agree closely with the other 
data for each Mach number. 

Shock-layer thickness for these conditions is presented in figure 11 for 
Mach numbers of 3 and 8.06. For each Mach number, the present calculated 
results agree closely with the other data. 

The shock-wave shape and sonic-line location about an ellipsoid with 
b/a = {2 are presented in figure 12 for Y = 1 -4 and M^ = 6.8. Experi- 
mental results from reference 11 are included for a comparison. The present 
computed results are in excellent agreement with the experimental results 
of reference 11. 


Hyperboloid 

The shock-wave shape and sonic-line location about a 39*37° hyperboloid 
at M^ = 6.8 and y = 1 *4 are compared with experimental data from refer- 
ence 11 (fig. 13). The present calculated results are in very good agreement 
with the experimental data. 
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Hemisphere Cylinder 


Pressure distributions on a hemisphere cylinder for y = 1-4 are presented 
in figure 14 for free-stream Mach numbers of 3, 4.03/ 5.06 f and 6.03. Experi- 
mental results from reference 14 are presented for comparison. The calculated 
results are in excellent agreement with the experimental data except at the 
highest Mach number (M^ = 6.03) where the results tend to deviate slightly at 
the downstream locations. The reason for this deviation is not known; however/ 
a similar trend was noted in reference 14. 

Shock-wave shapes for these conditions are presented in figure 15 where 
excellent agreement with the experimental data is shown. 


Spherically Blunted Cone 

The calculated shock-wave shape about a spherically blunted 30° half-angle 
cone in helium (y = 5/3) at M^ = 20.3 is presented in figure 16 and is com- 
pared with experimental data from reference 15. The agreement between the cal- 
culated and experimental data is excellent. The calculated drag coefficient for 
this case is 0.5589, which compares favorably with the measured value of 0.5540. 

The next case considered is a spherically blunted, 60° half-angle cone in 
air (y = 1.4) at a free-stream Mach number of 4.63. For this large half-angle 
cone, sonic conditions exist on the cone surface at the downstream corner. 

These conditions were imposed simply by fixing the corner surface pressure at 
the sonic value and carrying out the remainder of the calculation in the usual 
manner. The results for this case are presented in figures 17 and 18, along 
with experimental data from reference 16 for comparison. As illustrated in 
figure 17, the sonic line for this case extends beyond the last grid line so 
that the entire subsonic region is not enclosed within the computational space. 
Although this can lead to computational difficulties, it did not in the present 
case, and the results obtained are in good agreement with the experimental data. 
(See figs. 17 and 18.) 


Reverse Curvature Body 

The final case considered is a reverse curvature body with a nose shape 
(defined by eq. (22) ) matched to a 30° conical afterbody. This is a body shape 
that might result from ablation in the stagnation region due to high heating 
rates (as, for example, on a planetary probe). The converged coordinate system 
associated with this configuration is shown in figure 19; the shock-wave shape 
and sonic-line location are shown in figure 20; and the surface pressure distri- 
bution is shown in figure 21 . Since no known experimental data or theoretical 
results are available on this type body shape for comparison, these results are 
included here only to illustrate the type of complex body shape for which a flow 
field can be computed with the numerically generated coordinate system of refer- 
ence 6. In fact, it is with this type of complex body shape that it is antici- 
pated that the numerically generated orthogonal coordinate of reference 6 will 
be most useful. 
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It should be noted that a nonuniform spacing of grid points on the surface 
has been used in order to concentrate points in the expansion region near the 
beginning of the conical afterbody. This is a very useful feature of this coor- 
dinate system since any nonuniformity in the physical spacing of grid points is 
accounted for through the metric coefficients. 


CONCLUDING REMARKS 

A numerically generated orthogonal coordinate system (with the body sur- 
face and shock wave as opposite boundaries) has been applied with a time asymp- 
totic method to obtain steady-flow solutions for axisymmetric inviscid flow 
over several blunt bodies including spheres, paraboloids, ellipsoids, hyper- 
boloids, hemisphere cylinders, spherically blunted cones, and a body with a 
concavity in the stagnation region. Comparisons with experimental data and 
the results of other computational methods have demonstrated that accurate solu- 
tions can be obtained using this approach. The numerically generated orthogonal 
coordinate system used in the present paper should prove useful for applications 
to complex body shapes, particularly those with concave regions. In addition, 
the use of the present orthogonal coordinate system simplifies the form of the 
governing equations and simplifies the application of boundary conditions at the 
body surface and shock wave. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
January 22, 1980 
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APPENDIX 


CALCULATION OF PROPERTIES AT SHOCK WAVE 

A method similar to that presented in reference 9 has been used to calcu- 
late the shock velocity and other thermodynamic properties at the shock wave 
from the pressure. Consider the shock wave illustrated in sketch (a) which can 
be expressed as a function of Z,r s = r s (Z) 


Shock wave 


Sketch (a) 


The shock-wave angle $ s is given by the equation 


B s = tan -1 (dr s /dZ) 


and the components of free-stream velocity tangent and normal to the shock 
wave are given, respectively, by the equations 


(V t )eo = COS B £ 


( v n>co = v » sin 



(A!) 


(A2) 

(A3) 
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APPENDIX 


Changes in properties across the shock wave can be related to the normal 
component of free-stream velocity and the shock velocity, as illustrated in 
sketch (b) . 



Sketch (b) 


Conditions across the shock wave can be related by the following normal shock - 


wave equation (since = 

1) : 



(VnJeo - V s = p s [ 

<V n ) s “ 

V s] 

(A4) 

Poo + [< v n)oo - V s] 

1 2 

= Ps 

+ Ps [j v n) s " V S ] 

(A5) 

[< v n) OT - v s ; 

L\ Ia 

[(V n ) S - Vg] 2 

(A6) 

Ia>+ 2 


2 



Now, rearranging equation (A5) , the following expression can be obtained: 



(A7) 
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APPENDIX 


Combining this equation with equation (A4) and solving for - V S 1 , 

the following result is obtained 



Ps Poo 



(A8) 


Similarly equation (A6) becomes 


Is - 





(A9) 


Equations (A8) 



and (A9) can be combined with the ideal gas equation of state 
to yield the following result: 


1 

~s 


1 

+ "(PS * Poo) 



1 

" “(PS " Poo) 


(A10) 


Using this result, equation (A8) can be solved for the shock velocity V s 


V S - (V n )co - 



(All) 


Using equation (A4) , 


(V n ) s is obtained 


(V n )s = “ [(V n )co “ V s ] + V £ 


(A12) 
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APPENDIX 


Thus , the components of velocity downstream of the shock wave are simply 


( V t>oo - V oo COS 

(A13) 

(V n )s 

(AT 4) 


since the shock wave is the outer boundary of an orthogonal coordinate system. 
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Figure 3,- Orthogonal coordinate system in computational plane. 
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(a) M ffl = 6.8. 


Figure 6.- Shock-wave shape and sonic-line location for sphere, y = 1.4. 
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Figure 8.- Shock-layer thickness for paraboloid. 
Zp/Rjj = 1 .0; Y = 1.4. 
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Figure 9.- Pressure distribution on ellipsoid, b/a = 



(b) = 6.05. 

Figure 9.- Continued. 
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Figure 12.- Shock-wave shape and sonic-line location for ellipsoid. 

b/a = M= 6.8? y = 1 .4. 




Figure 13.- Shock-wave shape and sonic-line location for 39.37° hyperboloid 
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Figure 14.- Continued. 
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Figure 19.- Coordinate system on reverse curvature body. 
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